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ABSTRACT 

We propose a new method of gravitational wave de tection using a modi fied form 



of higher criticism, a statistical technique introduced bv lDonoho &: JinI (j2004l ). Higher 
criticism is designed to detect a group of sparse, weak sources, none of which are strong 
enough to be reliably estimated or detected individually. We apply higher criticism as 
a second-pass method to synthetic J-"-statistic and C-statistic data for a monochromatic 
periodic source in a binary system and quantify the improvement relative to the first- 
pass methods. We find that higher criticism on C-statistic data is more sensitive by ~ 6% 
than the C-statistic alone under optimal conditions (i.e. binary orbit known exactly) and 
the relative advantage increases as the error in the orbital parameters increases. Higher 
criticism is robust even when the source is not monochromatic (e.g. phase wandering 
in an accreting system). Applying higher criticism to a phase- wandering source over 
multiple time intervals gives a > 30% increase in detectability with few assumptions 
about the frequency evolution. By contrast, in all-sky searches for unknown periodic 
sources, which are dominated by the brightest source, second-pass higher criticism does 
not provide any benefits over a first pass search. 

Subject headings: gravitational waves — methods: data analysis — methods: statistical 
— pulsars: general — stars: binaries — stars: neutron 



1. Introduction 



Direct detection of gravitational waves appears likely in the near future. Existing terrestrial 
long-baseline interferometers, such as the Laser Interferometer Gravitational- Wave Observatory 
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(LIGO) and Virgo, have achieved their design sensitivity (jAbbott et alJl2009dl : lAccadia et al.ll2012l ) 
Next-generation interferometers now under construction are expected to detect tens of events per 
year, if contemporary estimates of compact binary coalescence rates are correct ([Abadie et al 



2010b|) 



Searches for periodic sources have the advantage of integrating over long observation times to 
increase the signal-to-noise ratio. The most likely periodic sourc es detectable by terrestrial inter- 



feroineters are rapidly rotatin g, slightly deformed neutron stars (jBildstenlll998l : lUshomirskv et al 



2000 



Melatos &: Payne |2005| ). They emit at the spin frequency and its first harmonic 2f^, 



(jjaranowski et al.l Il998ll. X-ray timing mea surements find 0.3 kHz < /* < 0.6 kHz for low- 



mass X-ray binaries (jChakrabartv et al.ll2003l ). placing these so urces squarely in the LIGO- Virgo 

j' 



2007al .l2008a. 2009a|, 



band. Directed searches fo r electromagnetically observed targets (jAbbott et al 
Abadie et aP boiOal. boil 



2004 



2007b. 



2008b': 



and blind, all-sky searches for unknown sources (jAbbott et al.l l2005i . 
have both been reported. 



Targeted searches for known pulsars use the radio ephemeris to guide the search and reduce 
computational expense by assuming the electromagnetic and gravitational wave phases track each 
other closely. One search for the Crab pulsar allowed for a mismatch of up to one part in 10 ^ 



between twice the radio pulse frequency and the gravitational wave frequency (jAbbott et al.ll2008bl ) . 
Targeted searches for the Crab and Vela pulsars have set direct upper limits on the wave strain 



(Abbott et al. 


2008b: 


' 

Abadie et al. 


2011 



A search of 78 pulsars using data 
from the third and fourth science runs of the LIGO and GEO 600 detectors set upper limits on the 
wave strain h and ellipticity e, the t ightest of which are h < 2.6 x 10^^^ for PSR J1603-7202 and 
e < 7.1 X lO-'^ for PSR J2124-3358 (jAbbott et al.ll2007bl ^. 



Blind, all-sky searches address more sources than targeted searches but are expensive computa- 
tionally as they cover a larger parameter domain to keep track of the unknown frequency evolution . 



A number of LIGO a ll-sky searches for periodic sources have been conducted (jAbbott et al 



2007a 



stein@Home project (jAbbott et al 



2005 



2008al . l2009a|HQL 

some of which leverage the distributed computing power of the Ein- 

3). Frequently, they ar e based on a maximum-likelihood 



2009a 



detection statistic, called the J^-statistic (jjaranowski et al.lll998l ). The J^-statistic is computed by 
coherently integrating over the ob servation time Tn h. t , assu ming a biaxial neutron star at a specific 
spin frequency and sky position. ICutler &: Schutzi (j2005l ) generalized the J^-statistic to apply to 
multiple detectors and sources. 

The computational expense of a fully coherent search for unknown sources becomes prohibitive 
as Tohs increases. Semi-coherent methods address this problem, dividing the observation time into 
intervals, which are individually searched coherently but combined incoherently. Semi-coherent 



(Wette 


2012 


)• 


Abbott et al. 


( 2008a ^ 



all-sky search for periodic sources. They described and compared three semi-coherent methods: 
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StackSlide, PowerFlux, and Hough num ber count (jBradv &: Creightonll200d: iKrishnan et al.ll2004l : 
Cutler et al.ll2005l : I Abbott et al.ll2008al ). For binary sources. [Messenger WoanI (|2007l ) derived the 
C-statistic, which incoherently combines individual J^-statistic values with a comb template that 
matches the orbital sideband structure. 

In this paper, we propose a new method for enhancing gravitational wave searches using a 
second-pass method known as h igher criticism. High er criticism was suggested originally by John 
Tukey and developed further bv lDonoho &: JinI (j2004l ). who proved that it has some mathematical 
optimality properties in the case of a specific noise distribution. The ultimate goal of any gravita- 
tional wave search is to identify an individual object as a definite source. Higher criticism is not 
suited for this purpose. It is designed to search against a known background for a group of sparse 
signals too weak to be detected individually. Higher criticism detects the presence of the group but 
cannot reliably estimate the waveform. We d iscuss applications of hi gher criticism to gravitational 
wave detection, focusing on the J-'-statistic ( Jaranowski et al.lll998l ) study. In general, 

we find that higher criticism is significantly more robust than other methods; it ignores some in- 
formation that other detection statistics incorporate, so may be less sensitive, but makes fewer 
assumptions about the form of the signal, especially important where the signal evolves during the 
search. 

The paper is structured as follows. In Section [21 we define the higher criticism statistic, 
compare the detectability of signals with different amplitudes and sparsities, and compute detec- 
tion thresholds. In Section [3l we briefiy review the J-'-statistic. In Section [U we apply a form 
of higher criticism to targeted searches for a binary pulsar and compare its performance with 
the C-statistic. We discuss how higher criticism can accommodate phase wandering in Section [5j 
While we apply higher criticism specifically to periodic gravitational wave searches in this paper, 
it is a general method with other application s, including detect ing non-Gaussianity in Wilkinson 
Microwave Anisotropy Probe (WMAP) data (jCavon et al.ll2005l ). 



2. Higher Criticism 



2.1. Theory 



The higher criticism statistic was introduced by iDonoho &: JinI (120041 ) to test the hypothesis 
that n independent and identically distributed (i.i.d.) samples Xi, . . . , Xn come from the same 
zero-mean Gaussian distribution A^(0, 1), against the alternative that a small fraction of them have 
a nonzero mean /U > 0, that is, to test: 

i.i.d. 



Ho 
Hi 



Xi ■ iV(0, 1 

i.i.d 



X., 



for z = 1, . . . , n , 

(l-e)iV(0,l) +eN{fj,,l) for i 



1, . . . , n . 



(1) 
(2) 



For example, Hq could correspond to a situation where all samples are just background noise with 
no signal (/x = 0), whereas in Hi a fraction e of samples contain signal {/j, ^ 0). Higher criticism is 
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particularly designed for the very challenging situation where the signal is sparse (e <^ 1) and weak 
{fj, is so small that the largest extreme values under Hi are essentially the same as those under Hq). 
Although we introduce the statistic in the context of Gaussian noise, it can be generalized easily to 
other distributions (see Section [2.3D. The case of correlated, or non- white, noise has been treated 



bv lHall k JinI (|20im 



Let Z denote a random variable with distribution A^(0, 1). The higher criticism statistic HC 
is computed from the n p-values pi = P[Z > Xi] generated by the n tests H^ i : Xi ^'^^ A^(0, 1) 
against Hi i : Xi N{n, 1), > 0, i = 1, . . . , n. It is defined by ( Donoho &: Jin 2004 ). 



HC = max 

l<i<n 



Vn[i/n-p(i)\ 



(3) 



where < p(2) ^ ■ ■ ■ ^ P{n) &re the piS sorted into increasing order. The statistic HC re- 
jects Hq when HC > g{n,a), where g{n,a) is a threshold chosen so that Pho[HC > g{n,a)] = 
-Pffn (reject Hq) < a, where Phq denotes the probability when Hq is true. [Our notation differs from 
Donoho &: JinI (j2004l ) in that we use the symbol g for the threshold instead of h to avoid confusion 



with the gravit ational wave strain h. 
HC > g{n, a 



asymptotically as n 
the threshold g{n, a) 



In other words, the higher criticism test detects a signal when 
Donoho Sz JinI (|2004l ) showed that, when Hq is true, one has HC ~ y^2 log log(n) 
oo for all a > 0, and that if < a < 1 and n is large enough, one can use 
i Y^2 log log(n). However for more general a and n, this does not work well 
in practice, as we illustrate in Section [2.4| where we discuss how to choose g{n,a) accurately as a 
function of n and a for finite n. 

Next, we quantify more precisely how small e and fi can be for the tes t based on HC to be 
able to distinguish Hq from Hi. Let r and /? be two positive parameters. iDonoho &: JinI (|2004l ) 
study the properties of the HC statistic for /u and e of the form 



e 



2r log n 

/3 



n 



(4) 
(5) 



These represent very difficult situations since, with r small, of this form is smaller than the mean 
of the upper extreme statistics, and with f3 large, the proportion of samples with non zero mean is 
very close t o zero. Clearly , a det ection is not possible if r is very small and /3 is very large at the 



same time. 



Donoho &: JinI (120041 ) showed that, as n — )• oo, there exists a function p(/3) defined by 



1/2 



1/2 < /3 < 3/4 , 
3/4 < /3 < 1 , 



(6) 



such that Hq and Hi are distingui shable only if r > p(l3 ). If r < /?(/3), there does not exist a test 



that can distinguish Hq from Hi. iDonoho &: JinI (j2004l ) proved that, under some conditions, the 



HC statistic is able to distinguish Hq from Hi throughout the detectable region r > p{f3), and 
that it has full asymptotic power, that is, P/f^ (reject Hq) — )• 1 as n — t- oo (here Phi denotes the 
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conditional probability when Hi is true). Another test has also been considered in the literature, 
namely the Neyman-Pearson like lihood ratio test, but it is less attractive since it requires r and (3 
to be known ( Donoho J iJ^2004V whereas the test based on the higher criticism statistic does not. 
Interestingly, iPonoho fc: Jin (,2004. ) also show that, in addition to the detection boundary, there is 
a second region r > j3 called the estimable region, where Hq and Hi can be distinguished a nd th e 
mean /i can also be estimated consistently. Figured! which can be found in lDonoho fc JinI ()2004l ). 
illustrates the detectable and estimable regions on the r-/3 plane. 



2.2. Distribution of HC 

In this section, we present briefly some results from Monte-Carlo simulations to illustrate the 
behavior of HC as a function of signal amplitude and sparsity e for n = 10^ sample values drawn 
randomly from the Hi distribution given by equation ([2]). 

Figure [2] displays a histogram of HC values for 10^ trials with e = 5 x 10"^ and // = 0, 0.1, 0.3, 1. 
Figure [3] displays a similar histogram with /i = 1 and e = 0, 5 X 10"^, 2 x 10"^, 5 x 10"^. The top left 
panel in both figures consists of samples drawn from the null distribution Hq. As /i or e increases, 
so too does the average value of HC. 



2.3. HC for a general noise distribution 



In Section 12.11 we defined HC in the context of Gaussian noise. As shown in Donoho &: Jin 



(120041 ) ■ higher criticism can be generalized to other noise and signal distributions, N and 5(/i) say, 
where describes the amplitude of the signal (e.g. mean of the Gaussian distribution in Section 
2.ip . In this more general setting we test to discriminate between the hypotheses 



Ho 
Hi 



i.i.d 



for 1 = 1, 



,n 



Xi ~ (1 - e)N + eS{fi) for i = 1, . 



(7) 
(8) 



In this context, the higher criticism statistic is computed as in equation but replacing the P(i)'s 
by the ordered values o f Pi = P\Z > Xj], i = 1, . . . ,n, where Z is a random variable that has the 



distribution Af; see also iDelaigle fc Halll (|2009l ). 

The detection boundary varies with the noise distribution. For example, when M is the 
distribution that we use later, the proportion e is defined as in equation ([5]), but the intensity of 
the signal is defined in terms of the noncentrality parameter of the distribution, through 



2r log n 



(9) 



instead of the mean fj. in equation ([!]). With these definitions of r and /3, the detection boundary 
is given by equation ([6]), as in the normal case. 
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Fig. 1. — Detectability regions on the r-(3 plane. Undetectable region, shaded gray; detectable 
region, where signal can be detected but not estimated, unshaded; estimable region, where signal 
can be both detected and estimated, shaded green. 
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Fig. 2. — Histograms of HC for 10^ Monte-Carlo trials, each constructed from n = 10^ samples 
drawn from distribution Hi with sparsity e = 5 x 10~^ and increasing amplitude /Lt. Upper left 
panel: noise only (fj, = 0), i.e. Hq distribution. Remaining panels: /x = 0.1,0.3,1 (upper right, 
lower left, lower right). 
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Fig. 3. — Histograms of HC for 10^ Monte-Carlo trials, each constructed from n = 10^ samples 
drawn from distribution Hi with amplitude /x = 1 and increasing sparsity e. Upper left panel: 
noise only (e = 0), i.e. Hq distribution. Remaining panels: e = 5 x 10~^,2 x 10~^,5 x 10"^ (upper 
right, lower left, lower right). 
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2.4. Detection thresholds 

In this section, we discuss the choice of the detection threshold g{n, a) for HC. As already 
noted earlier, the threshold satisfies Phq [HC > g{n, a)] = Phq (reject Hq) < a, and under somewhat 
restrictive conditions, it can be approximated asymptotically by y^2 log log(n). However, in finite 
samples and under more general conditions, this theoretical value y^2 log log(n) is not a good 
approximation to g{n,a). 

To illustrate this we conduct Monte Carlo simulations to determine the threshold g{n, a) for 
finite n. Note that, when Hq is true, the p- values are always independent and identically distributed 
according to the uniform distribution U{0, 1), and thus g{n,a) is independent of the specific noise 
distribution J\f. Therefore it suffices to generate p-values from the U{0, 1) distribution to compute 
the threshold. We run 10^ noise-only simulations with n = 10^, 10^, 10^, 10^ to determine Monte- 
Carlo thresholds. For each, we generate n p-values from U{0, 1) and compute HC. 

Figure H] displays the cumulative distribution function of HC for the simulated noise sam- 
ples. Threshold values obtained from the simulations by solving an empirical version of Phq{HC > 
g{n,a)) w a are listed in Tablefflfor n = 10^, 10^, 10^ 10^ and a = 0.5,0.1,0.05,0.01. For compar- 
ison, the theoretical value y^2 log log(n) is also quoted. From this and Figure HI we can see that 
this theoretical value generally does not work well in practice. In each case, a large fraction (> 
50%) of noise-only trials fall above this threshold. In the remainder of this paper, we commonly 
use the threshold g{n, 0.01) as determined from Monte-Carlo simulations with a = 0.01 (1% false 
alarm rate). 



3. J^-statistic 



The J^-statistic is an efficient detection statistic for periodic gravitational waves based on 
maximum likelihood. A num ber of blind and targeted searches for pe ri odic so ur ces have been carried 
out using it and its relatives (jAbbott et al.ll2004l . l2005l . l2007al . [2008al lbl. l2009al fl lAbadie et allboiOal . 



Table 1: HC threshold g{n,a). 

False alarm rate HC threshold g{n, a) 



a 


n = 10^ 


n = 10^ 


n = 10^ 


n = 10^ 


0.5 


2.10 


2.26 


2.37 


2.46 


0.1 


3.62 


3.66 


3.70 


3.73 


0.05 


4.72 


4.73 


4.73 


4.73 


0.01 


10.0 


10.1 


10.1 


10.0 


Theoretical threshold 


1.97 


2.11 


2.21 


2.29 
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Fig. 4. — Cumulative distribution function of HC when n = 10^ (red), n = 10^ (blue), n = 
10^ (green), and n = 10^ (black) samples drawn from Hq. Dashed vertical lines indicate the 
corresponding theoretical n — )• oo threshold g{n^ a) = y^2 log log(n). 
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201 ll ). We use the notation X^(A) to represent the noncentral (central if A = 0) distribution 
with k degrees of freedom and noncentrahty parameter A. In the absence of a signal, and assuming 
stationary, Gaussian noise, 2J^ is distributed according to a x'i(O) distribution, while in th e presence 
of a signal, 2J^ obeys a xlip'^) distribution (jjaranowski et al.lll998l : lAbbott et al.ll2007al ). In terms 
of the waveform, one writes 

where Sh{f) is the spectral noise density of the interferometer, Tobs is the total observation time, 
and h{t) = F+{t)h+{t) + Fx {t)hx (t) is the gravitational wave strain at the detector, which depends 
on the beam pattern functions and F^- For a biaxial star, the strains h^{t) and hx{t) depend 
on the wobble angle 9 between the rotation and symmetry axes, the inclination angle l between 
the rotation axis and line of sight, right ascension agkyi declination ^skyj polarization angle ^jJ, and 
the characteristic strain amplitude 



ho 



D 



(11) 



wher e e is the ellipticity, / is the moment of inertia, and D is the distance to the star (jjaranowski et al 
19981 ). We obta in an approximate relation between and fen by averaging equation ([TO]) over all 
relevant angles (jjaranowski et al.lll998l : IVigelius MelatosI 120091 ) . 



32 h?,T, 



Q^obs 



375 Sh 



(12) 



with the assumption that the spectral noise density is achromatic, viz. Shif) = Sh- Equation (jl2p 
allows us to translate detection limits on into limits on /ig throughout the rest of the paper. 

For a known source, the amplitude of the average s ignal that ca n be detected coherently with 
a 1% false alarm rate and a 10% false dismissal rate is (jWettell2012l ) 



ho = 15.6-1 



'^obs 



(13) 



Equation (jl3p is often seen with a factor of 11.4 instead of 15.6 (e.g., lAbbott et al.ll2007al ). The 
11.4 factor is obtained assuming the wobble angle 9 = vr/2, whereas we have averaged over 9 to 
arrive at equation (I13p . For an unknown source, we search across a large number of bins and Hq 
must be larger than the value in equation (jT3j) to rise above the background. 



Cutler fc Schutz (j2005l ) derived the J^-statistic for multiple detectors or multiple sources and 



found that combining the J-'-statistic values from multiple sources increases their detectability as 
long as the sources have a squared signal-to-noise ratio greater than one fifth of the brightest source. 
For example, of the known millisecond pulsars, it is advantageous to combine the five strongest 
ones but no more. 



Semi-coherent searches divide the observation time into intervals of length AT = Tobs/N. 
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The J-'-statistic is computed for each time interval i and combined to obtain 



N 



IT,, 



(14) 



As each 2Ti follows a xl distribution, 2Tsc follows a xIn distribution. 

Searches for periodic sources require detailed knowledge of the frequency evolution. For pulsars, 
the radio ephemeris is used to guide the search. Any difference between the radio and gravitational 
wave phases causes significant problems. For unknown sources, we ha ve even less information and 
are powerless to correct for phase wandering. Recently, ICutleii (|201ll ) proposed a "phase-relaxed" 
J-'-statistic for an all-sky search, modifying the semi-coherent detection statistic to accommodate 
a phase offset 6i for each coherent time interval i, while preserving a monochromatic phase model 
overall. 

In this paper we outline a new approach: applying higher criticism on a second pass to reanalyze 
the results of gravitational wave searches. Higher criticism relies on detailed knowledge of the 
background noise distribution. The detection statistics discussed so far assume stationary, Gaussian 
noise. Realistic detector noise is neither stationary nor Gaussian. We discuss this key point further 
in Section [6l For now, we persevere with the common, simplifying assumption of stationary, 
Gaussian detector noise in order to assess the overall viability of higher cr iticism, noting t hat 
modified forms of higher criticism have been cons tructed for correlate d noise (iHall &: JinI 120101 ) or 



when the noise distribution is imperfectly known (jPelaigle et al.l 12011 



4. Targeted binary search 



Low-mass X-ray binaries (LMXBs) are accreting neutron stars in a binary orbit with a low 
mass companion. X-ray pulsations and burs t oscillations place the sp in frequencies of discovered 
objects in the range 270 Hz < < 6 19 Hz ( Chakrabartv et al.l 120031 ). well below the centrifugal 
break-up frequency (jCook et al.lll994l ). Gravitational radiation is one effective means to balance 
the accretion torque and stall spin up (jBildstenlll998l ). In this section, we cor npare a second-pass 
search with higher criticism to the side-band C-statistic algorithm proposed by Messenger Woan 
(|2007l ). 



4.1. J^-statistic and C-statistic 



The gravitational wave signal from a neutron star in a b inary has its frequericy mo dulated into 
sidebands with spacing 1/P, where P is the orbital period ([Messenger WoanI 120071 ). Therefore, 
an J-'-statistic search sees a large number of relatively weak signals across many frequency bands. 
For a neutron star with intrinsic gravitational wave frequency /o, obeys a xiip^if)] distribution. 
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where the frequency dependent noncentrahty parameter can be written (jMessenger &: Woanll200?1 ) 



L^oJ 

PHf)=Pl E ^n(^o) W{f-fo,n: 
n=- [ZoJ 



(15) 



with Zq = 27r/oa, where [ZqJ represents the largest integer less than Zq, Jn is the n-th order Bessel 
function of the first kind, a is the light crossing time of the orbital radius projected onto the line of 
sight, /o,n = fo — n/P \s the frequency of the nth sideband around a source with frequency /o, Pq 
is the source signal-to-noise ratio squared [see equation (fTOj) ]. and W{f ) is the Fourier transform of 
the window function W{t), which equals 1 or when the detector is on or off respectively. Assuming 
no gaps in the data, one finds 



W{f) 



2 sin^(7r/rof,.) 



(16) 



obs 



With the window function defined as in [Messenger WoanI ( 20071 ). an extra factor of l/T^^^ 
included in equation (fT6|) so that C is dimensionless. 



IS 



The C-statistic sums incoherently the ^-"-statistic power in the orbital sidebands (jMessenger fc: Woan 



n/P) . 



(17) 



20071 ) as follows: 

C{f)= E 

n=-\_Zo\ 

For noise only, 2F is distributed according to a distribution. Hence, C, which is the sum of 
M 2T values, obeys a Y^ j^^ distribution, where M = 2[Zq\ + 1 is the total number of sidebands 
(jMessenger &: WoanI 120071 ). In the presence of a signal, 2F{f) is distributed according to a noncen- 
tral x\[p^{f)] distribution with a frequency dependent noncentra lity parameter p^f/) and C follows 
a noncentral X4jvf(-^) distribution with noncentrality parameter ( Messenger &: WoanI 120071 ) 



=-L^oJ 



n/P) 



(18) 



4.2. Higher criticism 

How does a C-statistic search, which assumes that P (and hence the sideband locations) are 
known, perform in comparison with a second-pass higher criticism search over J^-statistic or C- 
statistic values? We run Monte-Carlo simulations to answer this question. 

We simulate a targeted search for Sco X-1, similar to the one described in Sammut et al. (2013, 
in preparation) . The search parameters are P = 68023.84 s, a = 1.44 s, Tobs = 10 days, and source 
frequency /o = 400 Hz (assumed). We search over the range of frequencies 100 < / < 1000 Hz with 
frequency bin spacing 5f = l/(2Tobs) = 5.8 x 10^'' Hz, corresponding to n = 1.5 x 10^ frequency 
bins in total. The signal is modulated into M sidebands separated by 1/P. M depends on the 
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signal frequency, with 1811 < M < 18097 for 100 Hz < /o < 1000 Hz. The width in frequency 
space of the entire sideband structure (the comb) is 0.03 Hz < MjP < 0.27 Hz. 

Let us apply higher criticism to the problem. We take advantage of the relatively narrow 
sideband structure and divide the frequency domain [fmin, fmax] into windows of equal width 
w, each containing frequency bins, viz. [fmin, fmin + w]^ ifmin + w/2,frnin + 3w/2], [fmin + 
w, fmin + 2?/;], [fmax " Sw/2, fmax " w/2], [fmax " w, fmax]- The width w is choscu to be twice 
the maximum width of the sidebands, w = 2M/P, and windows overlap by 50% to ensure that the 
entire signal is contained entirely within a single window. A more sophisticated method would vary 
the width with frequency as the sideband width increases from 0.03 Hz for /o = 100 Hz to 0.27 Hz 
for /o = 1000 Hz, however, for simplicity, we construct our search windows based on the maximum 
sideband width at fmax = WOO Hz, obtaining w = 0.53 Hz, = 9.2 x 10^, and Ny^ = 3.4 x 10^. 

The data are synthesized as follows. First, the noncentrality parameter is computed for each 
frequency bin. The signal is modulated into sidebands and the frequency dependent noncentrality 
can be written (Sammut et al. 2013, in preparation) , 

[Zo] 2 
pHD^pI E JniZo)\w{f - fo,n)\ . (19) 
n=-lZo\ 

2 

We simplify this equation assuming W{f) 5{f) (true in the limit Tobs oo) and assume the 
J!^{Zq) all have similar amplitude, which we absorb into the constant Pq. Additionally, we adjust 
the frequency of each sideband /o,n to be equal to the corresponding closest frequency bin, nr(/o,n)) 
using nr(/) to represent the frequency bin closest to frequency /, neglecting the potential loss of 
amplitude caused by the mismatch between /o,n and nr(/o,n)- We denote the frequency on the i-th 
frequency bin using the 'prime' notation to differentiate between frequency bins and the exact 
frequency of the signal /q or the sidebands /o,n- The simplified noncentrality parameter used to 
generate the J^-statistic for each frequency bin f- is 

[Zoi 

P\fl)=pl E K[fl,nr{hn)] , (20) 

where K{x, y) = 5x,y is the standard Kronecker delta rebranded so that we can see the subscripts 
clearly. Equation (pO|) corresponds to p^{fl) equal to Pq for the frequency bin closest to each of 
the true sideband locations and zero elsewhere. We generate synthetic 2J^ values drawn from the 
x|[p^(/jO] distribution, for each frequency bin f^. 

To compute C-statistic values, we construct a comb centered at each frequency bin f^ assuming 
/o = fi, 

q{fl)= J2 K[finriflJ] , (21) 

n=-[ZU 
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with Z\ = 277 f -a and //„ = // + n/P. The C-statistic is computed by convolving the 2F{f) with 
q{f) (|Messenger k Woan|[2007l l 

C{f'^ = 2F*q{f^) . (22) 

We compute the higher criticism statistic as described in Section [2^ for the 2J^ and the C values 
contained in a given search window and call the results HC2T and HCc respectively. The noise 
distribution for 2T is xli^)- It is more complicated for C, as the number of sidebands = 2\_Z[\ +1 
depends on the frequency f[. Hence, the p-values for each C(/j-) are calculated with respect to the 
'X^Mi noise distribution (four degrees of freedom for each sideband summed). 

i 

To compare the performance of the three detection statistics C, HC2T and HCc, we compute 
thresholds with a false alarm rate a = 0.01. We have n 2T and C values, one for each frequency 
bin, and A'^ HC2T and HCc values, one for each search window. Thresholds are found by insisting 
that noise-only da ta fall below the threshold in every bin or window with probability at least 1 — a 



[e.g.- IWettd (1201J)], 

{P[C<CtH{a)]r>l-a, (23) 
{P [HC2T,C < 9{n^, a)]}''- > 1 - a , (24) 

where Cth{ot) is the C-statistic threshold. The HC threshold g{n^,a) is determined from Monte- 
Carlo simulations, as described in Section [2.41 We integrate the C-statistic probability distribution 
to find Cth{a) and use the result to compute the threshold noncentrality parameter Xth{a,S) for 
false dismissal rate 6 = 0.1 from 

/•oo 

/ dx Fix;4Ml^,^,0) = a/n , (25) 







dx F [x; 4M(^„,, Xth{a, 6)] = 6 , (26) 



where F{x;k,X) is the probability density function of the xlW distribution and AM[,^^^ is the 
maximum number of sidebands for all f'^ in the window. 

To simulate the gravitational wave search, we choose Pq and draw 2F values from the 'x\{p^{fi)\ 
distribution for each frequency bin f[. The search windows are constructed to ensure the entire 
sideband structure lies completely within a single window. Sidebands are also present in adjacent 
windows, however we focus only on a single window, assuming (conservatively) that all other win- 
dows contain only noise. Equation (j24l) makes the same assumption when computing the detection 
thresholds against which HC2j^ and HCc are compared. We conduct a 100 Monte-Carlo simula- 
tions for each and find the fraction where the computed statistic lies above its threshold value, 
which we call the detection rate. 

Figure [5] displays detection rates for C, HC2j^, and HCc as functions of increasing wave strain 
h/hth, where hth is defined as the wave strain corresponding to Xth- For 6 = 0.1 (detection rate 
of 90%), HCc detects wave strains 6% smaller than C alone and HC2j^ detects wave strains ~ 7 
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times greater than C. Additionally, higher criticism on 2J^ compared to 2J^ provides a advantage 
comparable to HCc over C (not displayed in Figure [5]). The results in Figure [5] assume perfect 
knowledge of the binary orbital period P. In summary, therefore, second-pass higher criticism 
on C slightly boosts the sensitivity of the C-statistic, while second-pass higher criticism on 2J^ is 
not competitive. Both higher criticism statistics, in particular HC2j^, are more robust than C, as 
demonstrated in the next section. 



4.3. Robustness of higher criticism 

One significant advantage of higher criticism is its robustness when applied to signals that differ 
from the assumed form, e.g. due to phase wandering. In this section, we explore the performance of 
the statistics C, HCc, and HC2T when the form of the signal remains unaltered but the true signal 
parameters differ from those assumed in the search. For a comb search targeting a binary source, 
the orbital period P determines the sideband spacing of the J^-statistic in the frequency domain. 
An error AP between the assumed and true value of P produces a comb that does not coincide 
with the actual sidebands. The C-statistic sums noise at each sideband whose frequency does not 
match a true sideband location, lowering the sensitivity. In contrast, HC2T only assumes that all 
sidebands are located within a particular frequency window. It does not rely on knowledge of the 
precise location of each sideband and consequently is unaffected by an error in P. Both HCq and 
the C-statistic lose sensitivity as AP increases, however one would expect the relative advantage of 
HCc compared to C is expected to remain. 



Observations of the orbital period of Sco X-1 have reported conflicting values f or P. iGottlieb et al 

1975 ) measured P = 68023. 8ib0. 09 s from archival optical observations. Recently. IVanderlinde et al 



2003) measured P = 68163.6 it 8.6 s with the Rossi X-ray Timing Explorer, but did not observe 



any signification periodicity at or near the previously reported value of 6 8023.8 s. Other LM XBs 



have larger uncertainty in their orbital parameters [e.g., see Tables 2-4 in lWatts et al.l ( 20081 )]. In 
light of the potential uncertainty in P, the robustness of higher criticism is valuable. 

We repeat the simulations described in Section [4.21 with the addition of an error AP between 
the assumed value of P used to construct the comb (?(/), and the true value used to generate the 
noncentrality parameters for the J^-statistic at every sideband. Figure [6] displays detection 

rates with C, HC2T-, and HCc as AP increases. As expected, HC2T is unaffected by AP because 
higher criticism does not use any information about the sideband locations. C and HCc both lose 
sensitivity as AP increases, but the relative advantage of HCc over C increases from 6% when 
AP = to 15% for AP = 1 s, 21% for AP = 3 s, and 29% for AP = 9 s. The difference 
between C and HC2T decreases from a factor of 7 when AP = to 5 for AP = 1 s, 3 for AP = 3 
s, and 2 for AP = 9 s. When AP = 0, the C-statistic peaks at the signal frequency /q. For 
AP > 0, the error between the assumed and actual location of the sidebands accumulates to the 
point where the outermost sidebands no longer coincide with the assumed comb. This reduces the 
total number of signal sidebands summed by the C-statistic, reducing its maximum amplitude (at 
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Fig. 5. — Monte-Carlo detection rate {1 — 6), equal to the fraction of simulations where the detection 
statistic lies above the detection threshold, as a function of wave strain, for the C-statistic (solid 
black), HCc (dashed red), and HC2T (dotted blue). 



- 18 - 




Fig. 6. — Monte-Carlo detection rate (1 — 6), equal to the fraction of simulations where the detection 
statistic lies above detection threshold, as a function of wave strain for the C-statistic (solid black), 
HCc (dashed red), and HC2T (dotted blue). The error AP between assumed and true binary 
orbital periods increases from top to bottom: AP = 0, 1, 3, 9 s. 
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/o). Additionally, the C-statistic peak broadens compared with AP = 0; this maximum amplitude 
requires less signal bins summed and can now be reached at a number of the central sidebands. This 
situation moves closer to multiple, equal-strength signals, the problem for which higher criticism is 
designed, explaining the increasing advantage of HCc over C. In contrast, the C-statistic requires 
a single source above threshold for detection. Interestingly, the sensitivity decrease of C and HCc 
stalls for 9 s < AP < 1000 s, when the accumulated error at the outermost sidebands is at least the 
sideband spacing. At this point, sidebands begin to overlap with other parts of the comb, shifted 
by integer multiples of the sideband spacing. As AP increases beyond 9 s, the number of sidebands 
correctly located by the template continues to decrease but is compensated for by a corresponding 
increase in the number of integer overlap sidebands. 



5. Phase wandering 

An accreting neutron star whose phase wanders in response to a variable accretion torque emits 
gravitational wave power in many frequency bins. Resampling is difficult as the phase model is 
usually unknown, e.g. there may be an offset between the radio/X-ray ephemeris and gravitational 
wave signal. Higher criticism can handle phase wandering robustly. 

In its simplest form, a semi-coherent J-"-statistic combines coherent time intervals of equal 
length according to, 

N 

J'sc{f) = Y.Hf) , (27) 

1=1 

where J-j is the usual .F-statistic for the i-th interval. This method assumes a monochromatic 



source or, with some modification, a source whose phase evolves in a known way. lAbbott et al 



(|2008al ) described three semi-coherent methods to search for periodic gravitational waves in LIGO 
data. The methods cannot be applied to a source whose phase wanders unpredictably, e.g., due 
to accretion torque. We apply higher criticism to reanalyze detection statistics from multiple time 
intervals with this situation in mind. 

We define nbins to be the number of frequency bins and ?igw to be the number of bins containing 
a gravitational wave signal, assuming for simplicity that all signal bins have equal amplitude. In 
practical terms, this situation corresponds to now distinct sources of similar strength, or new 
orbitally modulated sidebands of a source in a binary system. We also assume that nbins and new 
are the same in each time interval, therefore e = n-Gw/'^bins is constant. We allow for the signal to 
jump arbitrarily between bins (not necessarily adjacent ones) from one time interval to the next, 
e.g., due to phase wandering. 

For each interval, we obtain a 2J^ value for each of the ribins frequency bins. The ^gw signal 
bins follow a xlip'^) distribution, the remainder obey a xKO) distribution. Combining all 2J^ 
values for the intervals gives a total of n = A^^bins values. Converting from e and n to r and 
/3 through equations ([5]) and ([9]), we use the theoretical detection boundary [equation ^] to find 
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Phc(^^ '^GW) '^bins)) the minimum noncentrality detectable with HC over N time intervals, if racw 
of the nbins bins contain signal. This can be converted to a wave strain using equation (jl2p . The 
J-'-statistic threshold non centrality p% for a single interval is calculated similar to the C-statistic 



[equations ([25]) and ([26])] (jWettel 120121 ') 



dx F{x;4:,0) = ajn , (28) 



th 



dx F (x; 4, p%:) = <5i/"Gw _ (29) 



The semi-coherent J-'-statistic threshold pj-^^ can also be calculated from equations ()28p and (I29p . 
substituting 4A^, rather than 4, degrees of freedom and assuming a monochromatic source. The 
semi-coherent J^-statistic cannot be applied to source whose phase wanders unpredictably and the 
threshold reverts to pjr. 

Comparing detection thresholds for HC and J- is difficult because the detection boundary for 
HC, derived in the limit n — )• oo, is not equivalent to the J-'-statistic thresholds, computed for 
a specific false alarm and false dismissal probability. Furthermore, the HC detection boundary 
underestimates the threshold for finite n (see Section 12. 4p . We therefore recalibrate p\ju by a 
constant such that p%c(-^ ~ '^GW = !> '^bins) = Pjr, thereby arranging that both HC and T have 
equal power when = now = 1- To test this recalibration, we conduct Monte-Carlo simulations 
to determine pjj^ for a = 0.01, <5 = 0.1 when N = new = 1 for nbins = 10^, 10"^, 10^ 10^. We 
find that the Monte-Carlo p^^ agrees with pjr computed from equations ()28p and ()29p to within 
1%. However, it is important to mind the difference between each threshold when interpreting the 
following results. 

Figures [7] and [8] display required for detection with HC (thick curves) and J-" (thin solid 
and dashed curves) for different combinations of N, riQW) n-bins) and a = 0.01, 5 = 0.01. Figure 
[7] displays the thresholds pjj(^ and pjr as functions of new for (n-bins,-^) = (10^,1), (10^,20), 
(10^1), (10^20). Figure El displays as a function of N for (nbins, "Gw) = (10^1), (10^10), 
(10^, 1), (10^, 10^). For example, the search described in Section d] has nbins ~ 10^ and new ~ 10^ 
(the number of orbital sidebands). The semi-coherent J^sc is, in general, the most sensitive but 
as described above we are most interested in sources with unknown phase wandering, for which 
semi-coherent searches are not suitable. In all cases shown, HC outperforms J-. In Figure [71 we 
observe reduction in p'jj^^ by a factor of 1.5 to 2.9 for new = 10 and 2.3 to 7.0 for new = 100, 
compared to p'jjQ when new = 1. In Figure El we observe reduction in p'jj^j by a factor of 1.4 to 
2.2 for iV = 10 and between 2.1 and 3.0 for N = 100, compared to = 1. 

The robustness of HC described in Section [5] allows one to combine data from multiple time 
intervals to increased the sensitivity, even for a source whose phase wanders unpredictably. In 
principle, one might improve detectability while remaining robust by including some simple infor- 
mation about the source frequency evolution, e.g. limiting wandering to some physically motivated 
range. In Figures [7] and El we assume no correlation, allowing the unlikely possibility that the 




Fig. 7. — Noncentrality thresholds p\j(j (thick curves) and pjr (thin curves) as functions of new 
for (nbi„s,iV) = (10^,1), (106,20), (lO^l), (109,20) 
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frequency wanders across the entire range covered by the search. One simple method is to divide 
the frequency range into intervals wide enough to contain the source throughout the observation 
and calculate HC for each interval. 

6. Discussion 

Higher criticism is a recently formulated statistical method designed to detect the presence of 
a sparse collection of signals too weak to be detected individually. In this paper, we apply higher 
criticism as a second pass to reanalyze gravitational wave search statistics and explore the feasibility 
of applications in two contexts: a targeted binary search (e.g. LMXB) and a phase- wandering source 
(e.g., glitching pulsar). 

One advantage of higher criticism is its robust nature; it accommodates deviation from the 
expected phase evolution that hamper other search statistics. However, there is a trade-off. Higher 
criticism neglects some of the additional information used by other search methods which increase 
their performance under ideal conditions when the phase model is known. 

In Section HI we compare the performance of higher criticism and the C-statistic for a targeted 
binary search. The C-statistic is more sensitive than higher criticism applied to J^-statistic values by 
a factor ~ 7, however, higher criticism applied to C-statistic values gives a second-pass improvement 
of 6% over the C-statistic. Furthermore, higher criticism is more robust to an error AP between 
the true and assumed (from observation) binary orbital period. While the absolute sensitivity of 
both C and HCq decrease as AP increases, the performance of HCc relative to C increases to a 
sensitivity improvement of 16% for AP = 1 s and 29% for AP = 9 s. 

In Section [5l we consider a phase wandering source. The robustness of HC allows one to 
combine data from multiple time intervals and boost the sensitivity compared to a single interval, 
even for a unpredictable source. The noncentrality threshold decreases by a factor ^1.4 over 10 
intervals and a factor > 2 over 100 intervals. 

Another candidate for higher criticism is an all-sky search for unknown periodic sources. (We 
tested this in Appendix|A]). The results indicate that higher criticism provides no advantage over the 
J-"-statistic, the statistics perform similarly. This stems from the distribution of source amplitudes: 
if neutron stars are uniformly distributed across the Galaxy, the closest, strongest source dominates 
the detectability. There is no advantage using higher criticism, which is designed to detect a group 
of signals, when there is effectively only one signal present. 

While there is reason to be optimistic about possible applications of higher criticism, we draw 
attention to a number of significant concerns. Chief among these, higher criticism relies on knowing 
the background noise distribution. In this paper we adopted the common, simplifying assumption 
of stationary, Gaussian detector noise. In reality, detector noise is more complicated. However, 
it is well studied and modified forms of higher criticism have been developed for correlated noise 
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(jHall &: Jinll20ld ) or when the noise distribution is imperfectly known (jPelaigle et al.l 12011 



Additionally, higher criticism only detects the presence of a group of signals. It cannot directly 
identify an individual source, which remains the ultimate goal of gravitational wave detection. 
However, given the cheap computational cost of applying higher criticism as a second pass to 
reanalyze already computed search statistics, and the advantages related to its robustness, it has 
the potential to complement and enhance existing and future searches, especially as a guide to 
where to look harder in parameter space for just-too-weak sources. 
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A. All-sky search 

One might ask whether higher criticism can be applied pro fitably to all-sky s earches. After 



all, from an estimated Galactic population of 10^ neutron stars ( Arnett et al.lll989l ). only ~ 2000 



pulsars have been discovered as radio sources, of which 

current interferometers are most sensitive, /o > 100 Hz (iManchester et al 



1 0% have frequencies in t he range where 

Theories of 



neutron star quadrupoles suggest that only a handful of the observed pu lsars have any chance of 
being detectable with current sensitivity limits (jMastrano &: Melatosll2012l l. but many of the radio- 
quiet objects are much closer to the Earth than their radio-loud brethren. All-sky searches cover 
a range of frequencies at many sky locations, ho ping to detect one of these unknown so urces. A 
number of all-sky searches have been carried out (jAbbott et al.ll2005l . l2007al . l2008al . l2009al lbll3). but 
no detection has been announced. The searches typically use the J^-statistic, or related statistics. 
Hence they are candidates for reanalysis with higher criticism. 

In contrast to the other applications in this paper, the sources t argeted by an all - sky s earch 
have a power law (non-uniform) wave strain distribution. We follow ICutler &: Schutd (|2005l ) and 
assume neutron stars are spread uniformly throughout the Galaxy. Then the distance distribution 
separates into two parts: a local uniform three-dimensional distribution up to the thickness of the 
Galactic disk, and a uniform two-dimensional distribution beyond. Taking the disk to be 600 pc 
thick and 10 kpc in diameter, the 3D and 2D distributions correspond to D < 300 pc and 300 
pc < D < 5 kpc respectively, where D is the distance from Earth. We further assume that all 
sources have the same intrinsic amplitude, so that ho at Earth depends on distance alone. Changing: 
variab les from distance D to noncentrality parameter p^, [see equation (jlOp : also lCutler &: Schutz 
(| 20051 ) for details] we arrive at the distribution. 



(3D) , 
(2D) , 



(Al) 



where n2D and n^D are normalization constants which depend on ho, Shif) and Tabs through 
equation ([TO]) . 



To illustrate, we consider the sarn e search parameters as a completed LIGO all-sky search for 
periodic sources. lAbbott et al.l (|2007al ) searched the 10 hr of LIGO data with the best sensitivity 
from the second science run and computed the J-'-statistic for 1.6 x 10* frequency bins at each of 
3 X 10^ sky locations, a total of n = 5 x 10^^ J^-statistic values. We estimate the number of neutron 
stars with frequencies that fall in the range of the search (160-728.8 Hz) to be 10^ (of an estimated 
Galactic population of 10^). Under these assumptions, each of the 5 x 10^^ frequency bins is equally 
likely to contain some signal, with e = 2 x 10^^. 

We compute the Monte-Carlo higher criticism detection rate as a function of pfnaxl P%i where 
Pjr is the noncentrality threshold for a J-'-statistic detection with 1% false alarm rate and 10% 



^http://www.atnf. csiro.au/research/pulsar/psrcat 
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false dismissal rate [see equations ()28p and (j29p ]. We can then compare the noncentrality of the 
brightest source, p^axi required for detection with higher criticism with the equivalent quantity for 
a J^-statistic search. To save computation, we scale simulations from n = 5 X 10^2 to n' = 10^ by 
keeping parameters r and j3 constant in equations ([5]) and ([9]) and converting from e and for 
n = 5 X 10^^ to e' and p^' for n' = 10^. 

Signal values are drawn from a xKp^') distribution with noncentrality parameter p^' generated 
according to equation ()Aip . To fix n2D and nsD we specify p^^x fo^' ^ source at a distance of 10 
pc. The remaining F values, containing only noise, are drawn from the xKO) distribution, and 
we compute HC for the J^-statistic values. For each p^, we run 100 simulations and the detection 
rate equals the fraction of simulations with HC or T above their corresponding thresholds for 
detection. Finally, we convert noncentrality p^ to wave strain h through equation (jl2p to plot a 
more meaningful quantity. 

Figure [9] displays the Monte-Carlo detection rate as a function of hjhjr^ where h is the wave 
strain for a source at 10 pc, corresponding to p^axi normalized by hjr, the wave strain required for 
detection with the J^-statistic at 10 pc. Figure [9] shows that higher criticism provides no advantage 
over the J^-statistic; they perform almost identically. The fluctuations in Figure [9] relate to the 
randomly sampled source distribution. Despite the large number of sources, ct(p^) in equation (|Aip 
falls off too fast with p^ and effectively only the brightest source, rather than the collection of 
sources, is detected. To test this we repeated the above simulations using only the single strongest 
source (replacing the rest with noise) and found the same result. The similar performance of J- and 
HCjr agrees with Section [5l where Monte-Carlo simulations found HC and J- had equal thresholds 
when only one source was present. Therefore, for an all sky search, higher criticism does not provide 
any benefit over a first pass J-"-statistic search. 
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Fig. 9. — Monte-Carlo detection rate with higher criticism (thick dashed red) versus ^-"-statistic 
(thin sohd black) as a function of h/hj^, where h is the wave strain of a source located at 10 pc 
and hj: is the strain required by a source at 10 pc for detection with the J-'-statistic. 



